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Abstract 

This  paper  describes  extensions  of  goal-oriented  methods  for  a  posteriori  error 
estimation  and  control  of  numerical  approximation  to  a  class  of  highly-nonlinear 
problems  in  computational  solid  mechanics.  An  updated  Lagrangian  formulation 
formulation  of  the  dynamic,  large-deformation  response  of  structures  composed  of 
strain-rate-sensitive  elastomers  and  elastoplastic  materials  is  developed.  To  apply 
the  theory  of  goal-oriented  error  estimation,  a  backward-in-time  dual  formulation  of 
these  problems  is  derived,  and  residual  error  estimators  for  meaningful  quantities  of 
interest  are  established.  The  target  problem  class  is  that  of  axisymmetric  deforma¬ 
tions  of  layered  elastomer- reinforced  shells-of- revolution  subjected  to  shock  loading. 
Extensive  numerical  results  on  solutions  of  representative  problems  are  given.  It  is 
shown  that  extensions  of  the  theory  of  goal-oriented  error  estimation  can  be  devel¬ 
oped  and  applied  effectively  to  a  class  of  highly-nonlinear,  multi-physics  problems 
in  solid  and  structural  mechanics. 
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1  Introduction 


We  revisit  a  subject  addressed  nearly  a  half-century  ago  by  John  Argyris: 
’’Continua  and  Discontinua” ,  where  some  of  the  early  finite  element  approx¬ 
imations  of  a  large  class  of  nonlinear  problems  in  continuum  mechanics  were 
presented  [2],  Our  goal  here  is  to  bring  to  the  toolkit  a  new  methodology 
for  problems  in  nonlinear  continuum  mechanics:  goal-oriented  a  posteriori  er¬ 
ror  estimation  for  highly  nonlinear  dynamic  simulations  of  the  deformation  of 
submerged  bodies  subjected  to  shock  loading. 

Methods  for  developing  a  posteriori  estimates  of  finite  element  approxima¬ 
tions  of  linear  elliptic  boundary- value  problems  first  appeared  in  the  litera¬ 
ture  in  the  late  1970’s,  beginning  with  the  work  of  Ladeveze  [12]  for  elasticity 
problems  and  of  Babuska  and  Rheinboldt  [3]  on  two-point  boundary-value 
problems  and  followed  in  the  early  1980’s  by  extensions  to  elliptic  problems 
on  two-  and  three-dimensional  domains  [4,5,15].  Until  the  mid  1990’s,  except 
maybe  the  work  of  Gartland  [10],  virtually  all  of  the  methods  of  error  estima¬ 
tion,  an  essential  ingredient  in  mesh  adaptation  techniques,  were  applicable  to 
global  estimates  of  error  in  finite  element  approximations  of  linear  boundary- 
value  problems.  A  review  of  a  posteriori  error  estimation  can  be  found  in  the 
monograph  of  Ainsworth  and  Oden  [1],  Global  estimates  for  certain  classes 
of  nonlinear  elliptic  problems  were  contributed  by  Verfurth  [23]  see  also  [1], 
More  recently,  techniques  for  developing  a  posteriori  error  estimates  of  so- 
called  quantities-of-interest  which  are  functionals  of  solutions  of  linear  PDE’s, 
were  presented  by  Oden  and  Prudhomme  [16,17]  and  Becker  and  Rannacher 
(e.g.  [6]).  These  techniques  employ  optimal  control  strategies  and  involve  the 
solution  of  a  dual  problem  in  which  the  quantity  of  interest  appears  as  data. 
A  variety  of  applications  of  these  ideas  have  appeared  since  2000  (e.g.  [19— 
22]).  The  paper  of  Becker  and  Rannacher  [6]  appearing  in  2001  extended  the 
dual-based  theory  of  error  estimation  to  nonlinear  boundary-  and  initial- value 
problems,  and  the  paper  of  Oden  and  Prudhomme  [17]  extended  the  theory 
further  to  cover  estimation  of  both  modeling  and  approximation  error  in  2002. 
Further  extension  of  this  work  to  multi-scale  modeling  methods  is  described 
in  a  recent  report  [18]. 

While  the  developments  to  date  provide  an  abstract  mathematical  framework 
for  error  estimation  in  highly  nonlinear  problems,  no  applications  to  impor¬ 
tant  problems  in  nonlinear  continuum  mechanics  appear  to  have  been  made, 
owing  to  the  inherent  complexities  in  such  problems.  To  capture  features  of 
the  nonlinear  dynamics  of  solid  bodies  and  structures  under  shock  loading 
involves  a  host  of  complicated  features  and  has  been  the  focus  of  research  in 
computational  solid  mechanics  for  many  decades  (see  [14]  or  the  more  recent 
treatise  [7]).  The  analysis  of  the  evaluation  of  approximation  error  of  quan¬ 
tities  of  interest  in  such  applications  involves  solving  first  a  forward-in-time 
problem  for  the  system  response,  and  then  a  backward-in-time  problem  for 
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the  dual  solution  associated  with  the  particular  quantity  of  interest. 

In  the  present  investigation,  a  posteriori  error  estimates  for  key  quantities  of 
interest  are  derived  for  a  class  of  complex  and  highly-nonlinear  problems  in 
computational  solid  mechanics:  the  dynamical  behavior  of  a  heterogeneous, 
layered  shells  subjected  to  shock  loading.  The  models  considered  here  involve 
axisymmetric  deformations  of  thick  bodies-of-revolution  undergoing  very  high 
strains  and  strain  rates,  and  large  elastic  and  inelastic  deformations.  The  tar¬ 
get  applications  of  the  methods  developed  in  this  investigation  is  the  dynam¬ 
ical  behavior  of  elastomer-reinforced  steel  shells  subjected  to  high- intensity 
shock  loading. 

Following  the  introduction,  the  underlying  theory  of  a  posteriori  error  esti¬ 
mation  is  given  in  Section  2.  The  formulation  of  a  continuum  model  of  the 
problem  is  given  in  Section  3.  Here  a  framework  suitable  for  goal-oriented 
error  estimation  is  presented  and  algorithms  used  to  compute  the  solution 
of  the  discretized  space  and  time  problem  are  established.  Section  4  presents 
the  equations  of  goal-oriented  error  estimation  and  the  particular  solution 
technique  used.  The  method  of  assessing  the  fidelity  of  the  solution  is  also  dis¬ 
cussed.  The  geometry  and  data  of  the  target  application  problems  are  given  in 
Section  5  and  Section  6  presents  the  constitutive  equations  used  in  the  com¬ 
putational  model.  Detailed  numerical  results  are  presented  and  discussed  in 
Section  7.  Concluding  comments  are  collected  in  Section  8. 

2  Goal-Oriented  A  Posteriori  Error  Estimation  and  Control 

The  theory  of  goal-oriented  error  estimation  and  control  can  be  described  in 
terms  of  the  abstract  problem, 


Find  u  G  V  such  that 

B(u;v)  =  F(v)  Vv  €  V 


(1) 


where  £?(•,•)  is  a  semilinear  form,  nonlinear  in  the  first  entry,  v  is  a  test 
vector,  F(-)  a  linear  functional,  and  V  is  the  space  of  admissible  solutions, 
here  a  Banach  space  with  norm  ||  •  ||y.  Of  interest  is  the  value  of  a  functional 
Q  :  V  — *  M  at  solutions  a  to  (1);  the  quantity  of  interest.  The  problem  of 
determining  Q(u)  =  inf{Q(u)  :  v  E  V}  subject  to  (1)  is  an  optimal  control 
problem  characterized  by  the  pair  of  equations: 


Find  (u,p)  E  V  x  V  such  that 
B(u;v )  =  F(v)  Vv  €  V 
B'(u-,w,p)  =  Q(u;w)  \/w  G  V 


(2) 
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Here 


B\u;  w,p )  =  lim  -  [B(u  +  9w;p )  —  B(u;p)\ 
d^o  0 

Q\u-  w)  =  lim  ^  [Q{u  +  9w)  -  Q(u )] 

Problem  (2)i  is  the  dual  problem  associated  with  the  quantity  of  interest  Q 
and  the  primal  problem  (1)  (or  (2)2).  The  dual  problem  is  thus  linear  in  p, 
but  coupled  to  the  possibly  nonlinear  primal  problem. 

We  next  consider  a  family  {V^}  of  finite-dimensional  subspaces  of  V  with 
everywhere  dense  union 

U  vh 

/i-»  0 

generated,  for  example,  through  finite  element  approximations  of  functions 
(vectors)  in  V.  The  Galerkin  approximation  of  (2)  on  a  subspace  Vh  is  then 


Find  ( uh,ph )  eVhxVh  such  that 
B(uh;vh)  =  F(vh)  VvheVh 
B\uh ;  wh,  ph )  =  Q(uh;  wh)  Wwh  e  Vh 


(4) 


The  goal  is  to  estimate  the  approximation  error  £h  in  the  target  quantity  of 
interest  Q.  In  [17]  (see  also  [6]),  it  is  shown  that  to  within  terms  of  quadratic 
order  or  higher  in  the  error  components  eh  =  u  —  uh  and  £h  =  p  —  ph,  the  (a 
posteriori )  error  in  the  quantity  of  interest  is 

£h  =  Q(u)-Q(uh)^TZ(u\p)  (5) 

where  1Z(uh,  •)  is  the  residual  functional, 

n(uh-p)  =  F{p)-B(uh-p)  (6) 

We  note  that  the  Galerkin  approximation  ph  of  p  satisfies  the  orthogonality 
condition, 

n(uh-Ph)  =  o  (7) 

Various  goal-oriented  algorithms  may  be  constructed  for  systematically  reduc¬ 
ing  the  error  £h  [19-21],  Our  objective  is  to  formulate  the  field  equations  of 
nonlinear  continuum  mechanics  so  that  they  conform  to  the  structure  of  (2), 
and  to  then  develop  goal-oriented  methods  for  a  posteriori  error  estimation 
and  control  of  models  simulating  the  nonlinear  dynamics  of  layered  shell-like 
structures. 
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3  Variational  Formulation 


Basic  elements  of  the  formulation  and  much  of  our  notation  are  standard. 
We  use  an  updated  Lagrange  formulation  of  the  held  equations  of  nonlinear 
continuum  mechanics. 

3. 1  Weak  Form  of  the  Primal  Problem  The  forward  (primal)  problem 
for  an  updated  Lagrange  formulation  of  the  equations  governing  the  motion 
of  a  material  body  is  characterized  as  follows: 


Find  (u,  v)  e  V  such  that 

B((u,v);(z,w))  =  F((z,w))  V(z,w)eV 


(8) 


Here  V  =  Z  x  W  is  a  product  space  of  admissible  displacement- velocity  pairs 
and  £?(•;•)  and  F(-)  are  the  semilinear  and  linear  forms: 


B(( u,v);  (z,w))  = 

fT  r  dv 


f1  f  ,  drv  du  .  .  . 

/  /  (p—  ■  w  +  p—  •  z  —  pv  •  z  +  cr  :  Vxw)  dxdt 
J o  Jnt  dt  dt 

+  [  (p0v(X,  0)  •  w(X,  0)  +  p0u(X,  0)  •  z(X,  0))  dX  (9) 
J  r2n 


F((z,w))  = 

f  f  g  •  w  dAdi  +  /  [p0v0(X)  •  w(X,  0)  +  p0u0(X.)  •  z(X,  0)]  dX  (10) 
J o  J rf  JQo 

Here  we  ignore  body  forces  and  consider  the  motion  over  a  time  interval  [0,  T] 
of  a  material  body  occupying  a  current  configuration  ilt  C  M  at  time  t,  Q0 
being  the  reference  configuration.  The  energy  equation  is  not  considered  in  its 
weak  form.  In  (9)  and  (10), 

p  =  p(x,  t)  is  the  current  mass  density 
x  being  the  spatial  position 
of  material  particles  that  were  located 
at  position  X  in  the  reference  configuration 
v  =  v(x,  t)  the  velocity  held 
u  =  u(x,  t)  the  displacement  held 
cr  =  cr(x,t)  the  Cauchy  stress  tensor 

V^e,Athespatialgradient 

OXi 
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3.2  Time- Discretized  Formulation  The  governing  equations  with  ini¬ 
tial  conditions  for  the  Updated  Lagrangian  formulation  are  presented  below 


pj  =  Po 

p(X,0)  = 

Po(X) 

f  dv 

+  cr  :  V. 

/  Ph 7 ' w 

lnt  dt 

V(X,  0)  = 

v0(X) 

du 

~dt=Y 

u(X,0)  = 

u0(X) 

dA 


de  ^ 

pdt=a:D 
e(X,0)  =  eo(X) 


Here  the  strong  forms  of  the  mass  equation,  the  velocity-displacement  relation, 
and  conservation  of  energy  are  assumed.  Conservation  of  momentum  is  the 
only  equation  considered  in  its  weak  form. 


The  time  interval  [0,  T]  is  decomposed  into  subintervals  [fn,fn+1]  where  the 
time  step 

j.n+1  _  4 n  1 

At  = -  tn+2  =  ~(tn+1  +  tn) 

is  determined  by  the  Courant  condition.  The  algorithm  used  for  advancing 
the  velocity  and  displacement  fields  in  time  is  based  on  the  following  finite 
difference  scheme  for  the  acceleration 


dv 

dt 


t=tn 


yH,+  2  -  -yUl  2 


tn+  5  -  tn-\ 


A  leap  frog  method  is  used:  the  velocities  and  displacements  are  computed 
half  a  time  step  apart.  The  time-discrete  momentum  equation  is  then  of  the 
form 


/ 

J  eifn 


p"v"+5  .  w ndx 


=  / 


pnVn  2  .  w ndx 

[  gn  •  w ndA-  I  an  :  V,w" 

Jr%  Jn,n 


f  eifn 
+  (tn+i 


dx 


L  tn 


1) 


where  pn  is  calculated  from  the  mass  equation 


pnJn  =  po(X) 


The  displacement  is  updated  in  a  similar  manner  as 

un+i  =  un  +  ^n+1  _  jn^n+i 


(12) 
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While  the  energy  equation  is  approximated  point-wise  using  the  following 
finite  difference  scheme 


pnen+ 1  =  pnen  +  (i tn+1  -  tn)an  :  Dn"^  (13) 

and  is  only  used  to  update  the  yield  stress  of  the  materials.  As  the  name 
indicates,  during  the  Lagrangian  step  the  mesh  is  moved  with  the  material. 

Algorithm  1  explicit  time  algorithm  for  velocity  and  displacement 

use  initial  conditions  to  obtain  v5  and  u° 
for  n  =  1  to  ncycle  do 

determine  At  from  Courant  condition 
update  stress  crn  =  crn(un,  vn-2 }  en) 
use  (11)  to  compute  vn+2 
update  energy  using  (13) 
use  (12)  to  compute  un+1 
end  for 


The  specific  applications  of  this  formulation  to  be  considered  in  Section  5, 
6,  and  7  also  involve  contact  problems  encountered  in  layered  structures  in 
which  layers  are  not  bound  together  but  merely  are  in  physical  contact  in  the 
reference  configuration. 

The  layered  shells  are  therefore  modeled  with  a  sliding  interface  between  the 
layers,  with  the  interface  treated  as  a  contact  region.  The  standard  treatment 
of  contact  is  to  assume  that  the  bodies  in  contact  cannot  overlap  and  that  the 
tractions  of  the  surfaces  of  the  contact  region  satisfy  momentum  conservation 
at  the  interface.  In  the  computations,  these  two  conditions  impose  another  set 
of  nonlinear  equations  to  be  coupled  to  the  updated  Lagrangian  equations  of 
motion.  A  comprehensive  treatment  of  contact  and  implementation  may  be 
found  in  Belytschko  [7]. 

The  contact  algorithm  does  not  alter  the  dual  formulation.  From  the  view 
point  of  the  elastomeric  material,  in  which  the  dual  solution  was  computed, 
the  contact  region  is  merely  a  different  set  of  traction  boundary  conditions  on 
the  Neumann  boundary. 

4  The  Dual  Variational  Formulation 

There  are  two  distinct  sources  of  error  in  our  discretized  primal  problem.  One 
source  arises  from  the  finite  difference  approximation  of  the  time  derivatives 
and  the  second  source  of  error  from  the  finite  element  approximation  of  the 
solution.  This  analysis  is  primarily  concerned  with  the  latter.  In  the  following 
sections  our  problem  will  be  placed  in  a  suitable  framework  for  goal-oriented 
error  estimation  as  developed  by  Oden  and  Prudhomme  [17].  The  dual  problem 
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to  our  primal  problem  will  be  established  and  discretized.  Treatment  of  error 
estimation  is  restricted  entirely  to  the  elastomeric  layer  of  the  shell  structure. 
Consider  as  the  primal  problem  the  updated  Lagrangian  equations  formulated 
on  the  elastomer  layer,  see  Figure  2.  From  this  point  of  view,  the  remainder 
of  the  computational  simulation  will  merely  serve  to  prescribe  boundary  con¬ 
ditions. 

4-1  Dual  Problem  Procedures  for  deriving  the  dual  problem  has  been 
presented  in  Oden  and  Prudhomme  [16,17]  and  Becker  and  Rannacher  [6]. 
Consider  a  particular  quantity  of  interest  characterized  by  a  functional  Q, 
defined  on  V: 

<3((u,v))=  f  K(x)uzdx  (14) 

J  CIt 

where  K (x)  is  a  kernel  function  on  the  deformed  region  underneath  the  shock 
loading  material  in  the  final  configuration,  fR.  The  kernel  function  is  defined 
such  that  the  quantity  of  interest  represents  a  local  average  of  the  vertical 
displacement  uz  over  a  sub-domain  of  Hr  ■  The  relevance  of  this  quantity  of 
interest  will  become  clear  when  the  target  problem  is  presented  in  Section  5. 

Following  [17],  the  dual  problem  of  (8)  is  given  as  follows 


Find  (p,  q)  €  V  such  that 

p,q)) 


<2'((u,v);  (z,w))  V(z,  w)  G  V 


(15) 


where  B'  and  (}'  are  the  derivatives  of  B  and  Q.  respectively,  i.e. 

B'((u,v);  (z,w),  (p,  q))  = 

[5((u,v)  +  0(z,w);(p,q))  -R((u,V);(p,q))] 

and 

Q'((u,  v);  (z,  w))  =  lim  i  [Q((u,  v)  +  Q( z,  w))  -  Q((u,  v))] 

For  the  quantity  of  interest  particular  to  this  problem 

Q\( u,v);  (z,w))  =  /  K(x)zz(x)dx 

J  VL’J' 

Applying  the  definition  of  B’  to  our  semi-linear  form 
£'((u,v);  (z,w),  (p,  q))  = 

rT  r  d\V  dz 

/  /  (p-rr -q  +  vxq  :  Scru  :  Vxz  +  Vxp  :  Sav  :  Vxw-pw-p  +  p—  -p)  dxdt 

J o  Jnt  dt  dt 

+  /  (p0w(X,  0)  •  q(X,  0)  +  pDz(X,  0)  •  p(X,  0))  dX 


Where  Scrv  and  5cru  are  obtained  from  the  Taylor-expansions  of  the  Cauchy 
stress: 


Wj(Vx.u  +  9VX z,  Vxv  +  9\7X w)  =<rij(yx u,  Vxv)  +  9 


dcr 
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+  9— — -  w,„.„  +  (9(#2) 
<9(vm,n) 


with 


r  _n  _ 

ijmn 


da\j 

^(^,11) 


r  _ 

^  ijmn 


da\j 

d(vm,n) 


The  explicit  form  of  B'  shows  that  the  constitutive  equations  for  the  materials 
are  intimately  tied  to  the  error  estimate  calculations. 


Integrating  by  parts  to  transfer  the  time  derivatives  from  the  test  functions, 
(z,  w),  to  the  influence  functions,  (p,  q),  the  derivative  of  the  semi- linear  form 
becomes 


£'((u,v);  (z,w),  (p,  q))  = 

[  [  (-pw  •  ^  -  pz  •  ^  +  Vxq  :  8(7u  :  Vxz  +  Vxp 

Jo  Jnt  at  at 

+  [  (pw(x,  T)  •  q(x,  T)  +  pz(x,  T)  ■  p(x,  T))  dx 
JnT 


Scrv  :  Vxw  —  pw  •  p)  dxdt 


(16) 


From  (16),  it  is  seen  that  the  initial  conditions  of  the  dual  problem  involve 
initial  conditions  of  the  influence  functions  at  time  t  =  T.  This  means  that 
data  of  the  influence  functions  must  be  known  at  time  T  and  the  influence 
functions  are  propagated  backwards  in  time.  This  shows  that  any  algorithm  for 
computing  the  error  in  quantities  of  interest  requires  that  the  primal  problem 
is  first  solved  forward  in  time  from  0  to  T  and  the  dual  problem  from  T  to  0; 
i.e.  the  states  of  the  primal  problem  at  all  times  in  [0,  T]  must  be  available  to 
solve  the  dual  problem. 


4-2  Time  discretization  of  dual  problem  For  appropriately  chosen 
test  functions  the  dual  problem  may  be  written  as  two  coupled  first-order 
PDE’s  with  conditions  given  at  the  final  time: 


dq  ^ 
w  •  —  +  Vxp 

dt 


5<tv  :  Vxw  +  pw  •  p)  dx  =  0 


q(x,  T)  =  0 


/  (Vxq  :  5<tw  :  Vxz 

■hh 

Pi(x,T)  =  p2(x,  T)  =  0, 


Ps(x,T) 


)  dx  =  0 

p 


(17) 
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The  following  time-implicit  algorithm  is  implemented  to  reduce  the  amount  of 
data  needed  to  be  stored  while  managing  the  results  of  the  forward  problem 
in  coefficients  of  the  dual  problem  at  many  time  instances.  Note  that  the  time 
index  has  been  switched  from  ‘n’  to  ‘k’  to  emphasize  the  time  step  difference 
between  the  explicit  primal  problem  and  the  implicit  dual  problem. 


Derivatives  of  the  influence  functions  are  approximated  as 


dq 

dt 

dp 

dt 


1 


i  tk  —  tk  1 


t=r_2 


p k  -  pfc-x 

tk  -  tk~l 


Thus,  the  time  discretized  dual  problem  may  be  written  as 


;  (d<TU)fe-f  ;  VxZfe-t  -  pk~iz 


-o  •  V7  r,k-~ 


,Qk  i 

tk~  2 


—  dx  =  0  (18) 


tk  -  tk 

/ 

h  dx+ 


J  ^rr  +  pk~l •  pm 

tk~h  '  1 

[  Vxpk~ 2  :  (Scrv)k-^  :  dx  =  0 

I 


(19) 


where 


fc_i  pfc_1  +  pfc 

P  2  =  - - 


k  i  q  +  q 

q  2  =  » 


The  following  is  the  pseudo-code  used  to  implement  the  dual  problem. 

Algorithm  2  implicit  time  algorithm  for  influence  functions 
when  computing  primal  problem  save: 

(dcr“)fc-T  (5<rv)k~^,  zfc-T  wk~^ ,  and  tk  for  k  =  2, 3, ...,  Kf 
use  final  conditions  to  obtain  q Kf  and  pKf ,  where  tKf  =  T 
for  k  =  Kf  to  2  by  —1  do 

solve  linear  system  of  equations  resulting  from  (18)  and  (19) 
to  compute  qfc_1  and  pfc_1 

end  for 


4-3  A  Posteriori  Error  Estimation  and  Residual  Based  Refine¬ 
ment  Errors  in  the  quantity  of  interest  may  be  quantified  by  computation  of 
the  residual  [17].  The  residual  is  defined  as 

K((u '>'•);  (z,w))  =  F(( z,w))  -  B(( u‘,  v*);  (z,w)) 
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and  the  error  in  the  quantity  of  interest  is  approximated  by 


Q((u,v))-Q((uh,vh))^n((uh\vhl)-(p,q)) 

jhl,vAl);(ph2,qh2)) 

where  (rF1,  vh)  is  the  computed  solution.  For  these  calculations,  bilinear  shape 
functions  are  used  for  the  primal  solution,  (rF1 ,  vhl ) ,  and  second  order  serendip¬ 
ity  elements  are  used  for  the  dual  solution,  (p^2,  q^2).  Note  that  using  bilinear 
shape  functions  for  both  the  dual  and  primal  solution  would  result  in  a  residual 
of  zero  due  to  the  Galerkin  orthogonality  property. 


5  Application  to  Shock  Loaded  Layered  Shell 

The  geometry  of  the  problem  of  interest  is  given  in  Figures  2,3,  and  4.  A  shock¬ 
like  pressure  loading  is  applied  normal  to  the  entire  outside  surface  of  the  shell. 
The  pressure  loading  data  is  obtained  from  an  independent  simulation.  The 
time  evolution  of  the  pressure  loading  at  various  distances  from  the  centerline 
along  the  shell  is  given  in  Figure  1;  the  time  evolution  of  the  pressure  at  radii 
not  shown  may  be  taken  as  the  linear  interpolant. 


Figure  1.  Time  evolution  of  pressure  loading  on  shell  structure. 


Three  distinct  configurations  of  the  shell  are  considered.  In  the  first  configura¬ 
tion,  shown  in  Figure  2,  the  shell  is  composed  of  two  layers:  the  outer  layer  is 
steel  and  the  inner  layer  is  the  elastomeric  material.  As  indicated  in  Figure  1, 
the  pressure  load  is  applied  normal  to  the  entire  outer  surface  of  the  steel. 
The  intensity  of  the  loading  is  largest  within  a  radius  of  r  =  3.49cm  from  the 
centerline.  The  thickness  of  the  steel  is  tsteei  =  1.27cm  and  the  thickness  of  the 
elastomer  is  teiastomer  =  2.54cm.  Two  interface  conditions  between  the  steel 
and  elastomer  are  considered:  perfectly  bonded  and  frictionless  sliding. 
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In  the  second  configuration,  the  layers  are  reversed,  the  outer  layer  being  now 
made  of  elastomer  and  the  inner  layer  of  steel  (see  Figure  3).  The  pressure  load 
is  applied  to  the  outer  surface  of  the  elastomer.  Again,  the  thickness  of  the 
steel  is  tsteei  =  1.27cm  and  the  thickness  of  the  elastomer  is  teiastomer  =  2.54cm. 
Similar  to  the  first  configuration,  perfectly  bonded  and  frictionless  sliding 
interface  conditions  are  considered. 


Figure  2.  Geometry  of  the  problem  of  interest  for  the  first  configuration  of  the 
layered  shell. 


Figure  3.  Geometry  of  the  problem  of  interest  for  the  second  configuration  of  the 
layered  shell. 


The  final  configuration  considered  is  the  base  configuration,  the  elastomer 
is  removed  and  the  shell  is  composed  entirely  of  steel  (see  Figure  4).  The 
thickness  of  the  steel  is  tsteei  =  1.27cm.  This  configuration  will  be  used  as  the 
basis  for  comparison  of  the  effect  of  the  elastomeric  layer. 


pressure  loading 
steel 


Figure  4.  Base  configuration  of  shell. 
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All  geometries  are  assumed  to  have  a  rotational  axis  of  symmetry  as  illus¬ 
trated.  The  problems  are  modeled  with  a  Lagrangian  mesh  throughout  the 
domain.  The  updated  Lagrangian  equations  are  solved  throughout  the  mesh. 

Among  quantities  of  interest  are  those  that  indicate  the  effectiveness  of  the 
elastomeric  layer  as  a  reinforcement  of  the  steel  structure.  Therefore,  the  quan¬ 
tity  of  interest  in  this  investigation  has  been  defined  (14)  as  the  displacement 
of  a  region  of  the  layered  shell  under  the  shock  loading  source.  This  quantity 
will  be  used  as  a  measure  of  the  damage  inflicted  by  the  shock  loading. 

6  Constitutive  Equations 

The  Cauchy  stress  for  each  material  is  decomposed  into  its  dilatational  and 
deviatoric  components 


(T  —  —pi  +  S,  S,j  —  CTij  +  pSjj 

where  p  denotes  the  mechanical  pressure.  Thermodynamic  equilibrium  is  as¬ 
sumed  so  that  the  mechanical  pressure  equals  the  thermodynamic  pressure 
provided  by  an  equation  of  state. 


The  shock  conditions  suggest  that  an  equation  of  state  based  on  Hugoniot  data 
is  well  suited  for  this  problem.  Thus,  the  steel  is  modeled  as  a  hypoelastic- 
plastic  material  using  a  Mie-Gruneisen  equation  of  state  to  model  the  pressure. 


P  =Ph  +  poro(e-eh) 
PoC20r] 


Ph=Po  + 


(1  —  ST])2  ’ 


7]  =  1  — 


Po 

P 


eh  =  ea  +  Ph+Po 

^ Po 


(21) 


where  ph  is  the  Hugoniot  pressure,  e/,  is  the  Hugoniot  energy;  T0  is  the 
Gruneisen  gamma,  CQ  is  the  reference  sound  speed,  and  s  is  the  slope  of 
the  particle  velocity-shock  velocity  curve. 


The  deviatoric  stress,  S,  is  also  obtained  from  the  equation  of  state.  The  bulk 
modulus  is  proportional  to  the  partial  derivative  of  pressure  with  respect  to 
density  at  constant  internal  energy: 


K=p 


e 


The  Poisson  ratio,  u,  and  the  bulk  modulus,  K,  are  used  to  formulate  the 
deviatoric  component  of  the  stress: 


3K(1  -  2v) 
2(1  +  v) 


SVJ  =  GBdev 
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where  SVJ  is  the  Jaumman  rate  of  the  deviatoric  stress  tensor  and  G  is  the 
shear  modulus. 

Plasticity  is  modeled  following  classical  J 2  flow  theory  where  the  yield  stress, 
a y,  is  obtained  from  a  Johnson-Cook  plasticity  model.  All  material  data  pa¬ 
rameters  are  given  in  Tables  1  and  2. 

ov  =  (A+  Be ")  (1  +  Clnfe))  (1  -  17) 


Table  1.  Plasticity  parameters  for  steel 


A  (Pa) 

B 

n 

c 

m 

T*  (K) 

0.7922  x  109 

0.676 

0.26 

0.014 

1.03 

1793.0 

Here  e”  is  the  effective  plastic  strain,  ep  is  the  effective  plastic  strain  rate,  and 
T*  is  the  effective  temperature  defined  in  terms  of  the  melting  temperature, 
ambient  temperature,  and  Hugoniot  temperature;  Tm,  Tref,  and  7), ,  respec¬ 
tively. 


T*  = 


T  -T, 


ref 


T  -T 

±m  J-'T 


ref 


T  =  Th  + 


e~eh 


Based  on  the  experiments  performed  at  UCSD  [13],  under  shock  loading  con¬ 
ditions  the  deviatoric  component  of  stress  in  the  elastomer  is  insignificant 
compared  to  the  pressure.  Therefore,  the  elastomer  has  been  modeled  hydro- 
dynamically  for  these  simulations,  S  =  0.  A  Mie-Gruneisen  equation  of  state 
is  assumed,  equation  (21).  The  parameters  Ca  and  s  of  the  Mie-Gruneisen 
equation  of  state  are  obtained  from  a  graph  of  the  data  shown  in  Figure  5. 
The  rest  of  the  Mie-Gruneisen  parameters  are  given  in  Table  2, 

Table  2.  Mie-Gruneisen  EOS  parameters 


mat 

Co(f) 

S 

r0 

c^(kg  k) 

V 

steel 

7.85  x  103 

4.5  x  103 

1.49 

2.17 

4.48  x  102 

0.3 

elastomer 

1.07  x  103 

2.503  x  103 

2.05 

1.55 

1.173  x  103 

0.4998 

Assuming  a  hydrodynamic  state  of  stress  for  elastomer, 


Cij  P^ij 


and  neglecting  the  dependence  on  V.,,v.  the  stress  variation  in  the  dual  problem 
becomes 


S—V 

ucrijmn 


fsrru 

ijmn 


o 

doji 


dp  _  dp  po 

n/  \  ®ij  o  to  ^ mnO'i 

d{  um>n)  dp  J 2 
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Figure  5.  Elastomer  Data 
where  for  an  axisymmetric  problem 

C'll  =  (-fll  F22  M3, 3  -P33  —  F2i  F12  M3, 3  F33  +  Fn  F22  —  F 21  F12)  M2, 2 
+  -^11  +  -^11  M3, 3  F33 

C22  =  (Fi  1  F22  M3, 3  F33  —  F21  F\2  Uzt2,  F-,a  +  Fn  F22  —  ^21  -F12)  Mi,i 

+  F22  M3, 3  F33  +  F22 

C33  =  (-fll  M2, 2  F22  -F33  —  u2,2  ^21  ^12  ^33  +  ^11  Ftt)  Ml,! 

+  M2,l  F12  F33  +  M2, 2  F22  F33  +  My 2  -^21  M2,l  F\2  F33 
—  M2,i  Fn  Mi, 2  F22  F33  +  Mi, 2  F21  F33  +  F33 
Cl 2  =  (-F2I  F\2  M3, 3  F33  —  Fn  F22  M3, 3  F33  +  F2i  F12  —  -Fu  F22)  M2,l 
+  F 21  +  i*21  M3, 3  F33 

C21  =  (-^21  -^12  M3, 3  F33  —  Fn  F22  M3, 3  F33  +  F2 1  F12  —  Fn  F22)  mi,2 

+  F 12  M3, 3  F33  +  F 12 

(7i3  =  0  6*31  =  0  C23  =  0  C*32  =  0 


7  Simulation  Results 

Simulation  results  and  error  estimate  calculations  are  presented  in  this  section. 
Figure  6  illustrates  the  sequence  of  2-D  axisymmetric  quadrilateral  meshes 
used  to  model  the  problem  and  compute  error  estimates.  A  uniform  mesh 
was  used  for  the  elastomer  and  steel  directly  beneath  the  region  of  highest 
intensity  shock  loading;  the  mesh  was  grown  at  5%  away  from  this  region.  As 
a  measure  of  the  element  size,  N=4,  6,  7,  9,  11,  12  and  13  elements  were  used 
across  the  thickness  of  the  elastomer. 
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Figure  6.  Meshes 


7.1  Primal  Solution  Figure  7  presents  material  plots  of  a  simulation 
of  the  updated  Lagrangian  equations  of  motion  using  a  Lagrangian  mesh 
throughout  the  domain.  The  shell  is  composed  entirely  of  steel  and  had  no 
elastomeric  layers;  this  may  be  considered  the  base  model  from  which  com¬ 
parisons  can  be  made  between  configurations  that  contain  elastomeric  layers. 
Figure  8  presents  the  predicted  final  stress  state  of  the  shell.  Shown  is  de- 
viatoric  components  of  the  stress,  srr,  szz,  and  srz,  along  with  the  pressure. 
The  contours  are  given  in  units  of  Pascal’s  ranging  from  -650  MPa  to  +550 
MPa.  The  standard  sign  convention  is  employed:  tension  positive,  compression 
negative,  and  positive  traction  on  positive  element  face  is  positive  shear.  The 
plots  indicate  that  the  loading  is  not  localized  beneath  the  region  of  highest 
intensity  shock  loading.  Regions  away  from  the  centerline,  including  the  far 
edge  of  the  shell,  experience  large  stress  states.  The  contours  shown  are  in 
units  of  Pascal’s. 


Figure  7.  Results  for  Primal  Problem.  Base  model.  Steel  only,  top- left:  t  =  0  /is 
top-right:  t  =  100  ps  bottom-left:  t  =  300  ps  bottom-right:  t  =  500  ps 
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Figure  8.  Final  stress  states.  Base  model.  Steel  only,  top-left:  srr  top-right:  szz 
bottom-left:  srz  bottom-right:  pressure 


Figure  9  presents  material  plots  of  a  Lagrangian  simulation  with  the  config¬ 
uration  taken  as  steel  on  top,  elastomer  on  bottom,  and  a  sliding  interface 
between  the  steel  and  elastomer.  Figure  10  presents  the  predicted  final  stress 
state  of  the  shell  corresponding  to  Figure  9.  The  layers  of  the  shell  have  sepa¬ 
rated  near  the  centerline.  Because  of  the  separation  one  would  expect  a  final 
stress  state  similar  to  the  base  model.  However,  the  contact  region  near  the 
outer  edge  of  the  shell  appears  to  have  mollified  the  regions  of  high  shear 
stress  but  exacerbated  the  regions  of  high  radial  and  axial  stress  along  with 
the  pressure. 


Figure  9.  Steel  layer  on  top.  Sliding  Interface  top-left:  t  —  0  /is  top-right:  t  =  100 
/is  bottom-left:  t  =  300  /is  bottom-right:  t  =  500  /is 
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Figure  10.  Results  for  Primal  Problem,  final  stress  states  top-left:  srr  top-right:  szz 
bottom-left:  srz  bottom-right:  pressure 


Results  presented  in  Figure  11  are  for  the  same  configuration  as  those  in 
Figure  9  except  the  interface  between  the  steel  and  elastomer  is  modeled  as 
perfectly  bonded.  The  treatment  of  the  interface  condition  between  elastomer 
and  steel  is  not  seen  to  significantly  affect  the  final  deformed  shape  of  the  hull. 
However,  Figure  12  shows  that  the  final  stress  state  has  significantly  changed. 
The  loading  absorbed  by  the  steel  directly  under  the  region  of  highest  intensity 
shock  loading  appears  to  be  reduced  and  loading  away  from  the  centerline  has 
increased. 


Figure  11.  Steel  layer  on  top.  Bonded  interface  between  steel  and  elastomer,  top-left: 
t  =  0  ns  top-right:  t  =  100  fis  bottom-left:  t  =  300  fis  bottom-right:  t  =  500  ns 
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Figure  12.  Steel  Top.  Perfectly  Bonded,  final  stress  states  top-left:  srr  top-right:  szz 
bottom-left:  srz  bottom-right:  pressure 

Figure  13  presents  the  calculations  for  the  shell  configuration  with  the  elas¬ 
tomer  on  top  of  the  steel.  The  interface  between  the  steel  and  elastomeric  layer 
is  modeled  as  frictionless  sliding.  An  onset  of  a  pressure  instability  in  the  shell 
caused  by  the  contact  algorithm  combined  with  the  use  of  reduced  integration 
is  seen  at  time  t  —  100  /is.  The  stress  state  at  late  times,  Figure  14,  shows  the 
effect  of  this  instability.  Simulations  results  away  from  the  centerline  are  not 
as  reliable  as  the  results  near  the  centerline. 


Figure  13.  Elastomeric  layer  on  top.  Sliding  interface  between  steel  and  elastomer, 
top-left:  t  =  0  fis  top-right:  t  =  100  jis  bottom-left:  t  =  300  /is  bottom-right:  t  = 
500  jis 
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Figure  14.  Elastomeric  layer  on  top.  Sliding  interface  between  steel  and  elastomer, 
final  stress  states  top-left:  srr  top-right:  szz  bottom-left:  srz  bottom-right:  pressure 

The  results  from  the  layer  configuration  with  the  elastomer  on  top  and  the 
interface  between  the  shell  layers  modeled  as  perfectly  bonded  are  presented  in 
Figure  15.  Notice  that  the  final  stress  state,  Figure  16,  is  surprisingly  similar 
to  the  stress  state  shown  in  Figure  12;  indicating  that  the  interface  between 
the  layered  shell  is  more  significant  than  the  position  of  the  layers. 


Figure  15.  Elastomeric  layer  on  top.  Bonded  interface  between  steel  and  elastomer, 
top-left:  t  =  0  fis  top-right:  t  =  100  jis  bottom-left:  t  =  300  fis  bottom-right:  t  = 
500  fis 
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Figure  16.  Results  for  Primal  Problem,  final  stress  states  top-left:  srr  top-right:  szz 
bottom-left:  srz  bottom-right:  pressure 

1.2  Dual  Solution  The  dual  solution  corresponding  to  the  results  given 
in  Figure  9,  steel-on-top,  sliding  interface,  are  shown  in  Figures  17,  18,  and 
19.  The  specific  quantity  of  interest  is  the  displacement  of  a  region  of  the 
elastomeric  layer  under  the  loading  of  the  region  of  highest  shock  loading 
intensity.  Figures  17  and  18  show  the  radial  and  axial  components  of  the  dual 
solution.  Time  t  =  90//. s  of  the  axial  dual  solution  illustrates  the  position  of  the 
region  used  for  the  error  computation.  The  magnitude  of  the  axial  component 
of  the  dual  solution  is  seen  to  increase  as  the  integration  backwards  in  time 
proceeds. 


— Se-08 


E6 


Figure  17.  Axial  Component  of  Dual  Solution.  Steel  layer  on  top.  Sliding  interface, 
top-left:  t  =  90  fis  top-right:  t  =  70  fis  bottom-left:  t  =  40  fis  bottom-right:  t  =  0 

flS 
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Figure  18.  Radial  Component  of  Dual  Solution.  Steel  layer  on  top.  Sliding  interface, 
top-left:  t  =  90  fis  top-right:  t  =  70  jis  bottom-left:  t  =  40  fis  bottom-right:  t  =  0 
/is 

The  quantity  of  interest  was  chosen  as  the  average  displacement  in  the  axial 
direction,  consequently,  the  radial  component  of  the  dual  solution  is  seen  to 
be  an  order  of  magnitude  smaller  than  the  axial  component.  In  terms  of  the 
final  error  estimate  calculations,  this  means  that  forces  and  stresses  affecting 
the  radial  momentum  of  the  elastomer  have  a  significantly  less  contribution 
to  the  final  error  estimate. 
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Figure  19.  Spatial  Distribution  of  Residual.  Steel  layer  on  top.  Sliding  interface.  N= 
#  of  mesh  cells  across  thickness  of  elastomer,  top-left:  N  —  6,  top-right:  N  =  10, 
bottom-left:  N  —  14,  and  bottom-right:  N  =  18 

Figure  19  plots  the  spatial  distribution  of  the  contributions  associated  with 
each  node  over  the  entire  time  interval  to  the  error  in  the  quantity  of  interest 
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for  different  meshes.  The  region  where  the  magnitude  of  these  contributions 
is  largest  indicates  where  mesh  refinement  is  needed.  The  regions  close  to  the 
region  of  interest  and  near  the  highest  intensity  loading  influence  the  quantity 
of  interest  the  most.  Thus  mesh  refinement  only  directly  underneath  the  region 
of  highest  shock  loading  intensity  is  needed  for  accurate  computations  of  the 
quantity  of  interest. 
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Figure  20.  Axial  Component  of  Dual  Solution.  Steel  layer  on  top.  Bonded  interface, 
top-left:  t  =  90  /is  top-right:  t  =  70  fis  bottom-left:  t  =  40  fis  bottom-right:  t  =  0 
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Figure  21.  Radial  Component  of  Dual  Solution.  Steel  layer  on  top.  Bonded  interface, 
top-left:  t  =  90  /is  top-right:  t  =  70  /is  bottom-left:  t  =  40  /is  bottom-right:  t  =  0 
Us 

Figures  20,  21,  and  22  illustrate  the  dual  solution  results  corresponding  to  the 
steel-on-top  configuration  of  the  layered  shell  with  a  bonded  interface  between 
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the  layers.  Figures  20  and  21  are  the  axial  and  radial  components  of  the  dual 
solution.  The  results  are  similar  to  before.  The  regions  near  the  quantity  of 
interest  and  where  the  highest  intensity  loading  originates  require  the  most 
refinement.  The  majority  of  the  error  contributions  are  a  consequence  of  the 
loading. 
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Figure  22.  Spatial  Distribution  of  Residual.  Steel  layer  on  top.  Bonded  interface. 
N=  #  of  mesh  cells  across  thickness  of  elastomer,  top-left:  N  =  6  top-right:  N  =  10 
bottom-left:  N  =  14  bottom-right:  N  =  18 

Shown  in  Figure  23  is  a  graph  of  the  global  contributions  in  space  at  each  time 
step  to  the  error  in  the  quantity  of  interest  for  both  the  sliding  interface  and 
bonded  interface  conditions.  The  time  instances  that  contribute  the  most  error 
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Figure  23.  Temporal  Distribution  of  Residual.  Steel  layer  on  top. 

correspond  to  the  loading  of  the  elastomer  by  the  pressure  wave.  Figure  24 
plots  the  evolution  of  the  final  error  estimates  as  a  function  of  time. 
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Steel  on  Top 


Figure  24.  Temporal  Evolution  of  Final  Error  Estimate.  Steel  layer  on  top. 

Figures  25,  26,  and  27  are  the  dual  solution  results  corresponding  to  Figure  13, 
elastomer  on  top,  sliding  interface.  The  radial  and  axial  components  of  the  dual 
solution  are  shown  in  Figures  25  and  26.  Note,  for  this  case  the  dual  solution 
was  only  computed  for  T  =  15  ps.  Large  displacement  gradients  at  late  times 
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Figure  25.  Axial  Component  of  Dual  Solution.  Elastomeric  layer  on  top.  Sliding 
interface,  top-left:  t  =  13  jis  top-right:  t  =  10  /is  bottom-left:  t  =  7  /is  bottom-right: 
t  =  1  /is 

cause  numerical  instabilities  in  the  dual  solution  and  are  the  cause  of  the 
time  interval  restriction.  However,  the  results  are  similar  to  those  obtained 
before.  The  regions  in  the  mesh  where  the  highest  intensity  structural  loading 
originates  require  the  most  refinement  and  the  residual  is  seen  to  significantly 
increase  for  higher  mesh  resolutions. 
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Figure  26.  Radial  Component  of  Dual  Solution.  Elastomeric  layer  on  top.  Sliding 
interface,  top-left:  t  =  13  fis  top-right:  t  =  10  fis  bottom-left:  t  =  7  fis  bottom-right: 
t  =  1  fis 
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Figure  27.  Spatial  Distribution  of  Residual.  Elastomeric  layer  on  top.  Sliding  inter¬ 
face.  N=  #  of  mesh  cells  across  thickness  of  elastomer,  top-left:  IV  =  6  top-right: 
iV  =  10  bottom-left:  IV  =  14  bottom-right:  IV  =  18 


The  dual  solution  results  corresponding  to  Figure  15,  elastomer  on  top,  bonded 
interface,  are  shown  in  Figures  28,  29,  and  30.  The  radial  and  axial  compo¬ 
nents  of  the  dual  solution  are  shown  in  Figures  28  and  29.  Again,  the  dual 
solution  was  only  computed  for  T  —  15  /js.  The  results  are  consistent;  the  axial 
component  of  the  dual  solution  is  an  order  of  magnitude  larger  than  the  ra¬ 
dial  component  and  regions  in  the  mesh  where  the  highest  intensity  structural 
loading  originates  contribute  significantly  to  the  error. 
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Figure  28.  Axial  Component  of  Dual  Solution.  Elastomeric  layer  on  top.  Bonded 
interface,  top-left:  t  =  13  fis  top-right:  t  =  10  fis  bottom-left:  t  =  7  fis  bottom-right: 
t  =  1  fis 

-3e-10 

-2.4e-10 

-1.8e-10 

—  1.2e-10 

— 6e-11 

-  0 

— 6e-11 

— 1.2e-10 

-1.8e-10 

-2.4e-10 

-3e-10 

Figure  29.  Radial  Component  of  Dual  Solution.  Elastomeric  layer  on  top.  Bonded 
interface,  top-left:  t  =  13  fis  top-right:  t  =  10  fis  bottom-left:  t  =  7  fis  bottom-right: 
t  =  1  fis 
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Figure  30.  Spatial  Distribution  of  Residual.  Elastomeric  layer  on  top.  Bonded  inter¬ 
face.  N=  #  of  mesh  cells  across  thickness  of  elastomer,  top-left:  N  =  6  top-right: 
N  —  10  bottom-left:  N  —  14  bottom-right:  N  —  18 
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The  temporal  contributions  to  the  residuals  as  a  function  of  time  for  both 
sliding  and  bonded  interface  conditions  are  shown  in  Figures  31  and  32.  For 
this  short  time  interval  the  quantity  of  interest  is  seen  to  being  highly  effected 
by  the  loading  throughout  the  duration  of  the  simulation. 
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Figure  31.  Temporal  Distribution  of  Residual.  Elastomeric  layer  on  top. 


Figure  32.  Temporal  Evolution  of  Final  Error  Estimate.  Elastomeric  layer  on  top. 

Figures  33  and  34  show  final  error  estimates  corresponding  to  equation  (20). 
The  error  estimate  is  plotted  against  the  number  of  degrees  of  freedom  in 
the  mesh.  As  expected  the  error  estimates  are  seen  to  decrease  with  mesh 
resolution  and  the  error  estimates  for  the  configuration  with  elastomer  on  top 
is  seen  to  have  less  accumulation  of  error  because  of  the  shorter  time  interval 
that  the  problem  was  run  on. 

The  interface  condition  between  the  steel  and  elastomer  for  the  configuration 
with  steel  on  top  is  seen  to  produce  a  substantial  difference  in  the  rates  of 
convergence.  This  may  be  attributed  to  the  difference  in  the  loading  of  the 
elastomer.  For  the  configuration  with  elastomer  on  top,  the  rates  of  conver¬ 
gence  are  roughly  the  same  for  both  interface  conditions.  Over  the  short  time 
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interval  in  which  error  estimates  were  computed  the  initial  loading  for  both 
interface  conditions  is  approximately  the  same. 


Figure  33.  Estimate  Error  as  a  function  of  number  of  degrees  of  freedom.  Steel  layer 
on  top. 


Figure  34.  Estimated  Error  as  a  function  of  number  of  degrees  of  freedom.  Elas¬ 
tomeric  layer  on  top. 

8  Conclusions 

In  summary,  a  new  simulation  tool  has  been  developed  to  study  nonlinear  dy¬ 
namics  and  large  deformations  of  shock-loaded  structures.  The  theory  of  goal 
oriented  error  estimation  has  been  applied  to  a  class  of  highly  nonlinear  shock 
loaded  problems.  Material  properties  of  the  elastomer  were  exploited,  allowing 


29 


approximations  to  the  Cauchy  stress,  resulting  in  a  computationally  feasible 
dual  problem  for  highly  complex  phenomena.  The  dual  solution  was  used  to 
obtain  converging  error  estimates  for  a  meaningful  quantity  of  interest.  This 
work  may  be  used  as  a  foundation  for  developing  adaptive  meshing  algorithms 
for  this  class  of  nonlinear  problems. 

In  this  study,  extensive  numerical  calculations  were  performed  on  various  con¬ 
figurations  of  layered-shell  structures.  The  primary  effect  of  the  elastomer  was 
to  redistribute  the  regions  of  high  stress  concentration  through  out  the  steel. 
Results  also  showed  that  the  differences  in  the  final  stress  state,  as  compared 
to  the  base  model,  is  dominated  by  the  interface  between  the  layers  of  the  shell; 
i.e.  the  positions  of  the  layer  is  not  significant,  but  the  mechanism  that  allows 
the  transfer  of  energy  and  momentum  between  the  layers  is  very  important. 
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